home *** CD-ROM | disk | FTP | other *** search
/ Eagles Nest BBS 8 / Eagles_Nest_Mac_Collection_Disc_8.TOAST / Developer Tools⁄Additions / MacScheme20 / Mathlib / ratfun.scm < prev    next >
Encoding:
Text File  |  1989-04-27  |  5.7 KB  |  182 lines  |  [TEXT/????]
  1. ;;; $Header: ratfun.scm,v 1.6 88/09/21 13:38:15 GMT hal Exp $
  2. ;;;; 6.003 Rational Functions
  3.  
  4. (if-mit
  5.  (declare (usual-integrations = + - * /
  6.                  zero? 1+ -1+
  7.                  ;; truncate round floor ceiling
  8.                  sqrt exp log sin cos)))
  9.  
  10. ;;; This file contains these definitions:
  11. ;;;
  12. ;;; (procedure->rat-func <procedure> <degree>) -> ( numer denom )
  13. ;;;
  14. ;;; (rat-func->procedure <ratfunc> [procedure])
  15. ;;;
  16. ;;; (rat-func-compose <rf1> <rf2>)
  17.  
  18. ;;; The following procedure, given a function of complex frequency
  19. ;;; derived from a circuit of known degree, returns a list of lists of
  20. ;;; numbers that are the coefficients of the numerator and denominator
  21. ;;; polynomials.  The list of lists ((a0 a1 a2) (b0 b1 b2 b3))
  22. ;;; corresponds to the rational function
  23. ;;;
  24. ;;;                  a2*s^2 + a1*s + a0
  25. ;;;              --------------------------
  26. ;;;              b3s^3 + b2*s^2 + b1*s + b0
  27. ;;;
  28. ;;; The degree supplied, n, corresponds to the max degree of the
  29. ;;; denominator or (which is the same thing) to 1 + the max degree of
  30. ;;; the numerator.  If n is too small, the procedure will return the
  31. ;;; rational function that agrees with the values of the given
  32. ;;; procedure at 2n test points generated by TEST-POINTS.  If n is too
  33. ;;; large, an error is signalled, indicating that the equations are
  34. ;;; not independent.
  35. ;;; The answer is normalized so that the highest non-zero power of s
  36. ;;; in the denominator has coefficient 1. 
  37.  
  38. (define (procedure->rat-func procedure n)
  39.   (if (< n 1) (error "Bad maximum degree -- PROCEDURE->RAT-FUNC" n))
  40.   (let* ((coeffs (rat-fun-vector procedure n))
  41.      (a-part
  42.       (generate-list n
  43.          (lambda (i)
  44.            (heuristic-round-complex (vector-ref coeffs i)))))
  45.      (b-part
  46.       (generate-list n
  47.          (lambda (i)
  48.            (heuristic-round-complex (vector-ref coeffs (+ i n)))))))
  49.     (list (reverse (truncate-leading-zeros (reverse a-part)))
  50.       (append b-part '(1))))) ;b[n]=1
  51.  
  52.  
  53. ;;; Inverse of `procedure->rat-func' is easy:
  54.  
  55. (define (rat-func->procedure r . optional)
  56.   (let ((r (check-rat-func r)))
  57.     (let ((point (if (null? optional) identity-operator (car optional)))
  58.       (numerator (car r))
  59.       (denominator (cadr r)))
  60.       (/ (poly-value numerator point) (poly-value denominator point)))))
  61.  
  62. ;;; Below this point we implement the PROCEDURE->RAT-FUNC procedure.
  63.  
  64.  
  65. (define (truncate-leading-zeros l)
  66.   (cond ((null? l) (list 0))
  67.     ((= (car l) 0) (truncate-leading-zeros (cdr l)))
  68.     (else l)))
  69.  
  70.  
  71. ;;; We generate the 2n by 2n matrix M of coefficients of the a[i]
  72. ;;; and b[i] obtained by evaluating f at 2n points and setting
  73. ;;; this equal to the desired rational form.  
  74.  
  75. (define (rat-fun-vector f n)
  76.   (let ((e (assemble-matrix-equation f n)))
  77.     (let ((A (car e)) (B (cadr e)))    ;AX=B
  78.       (let* ((lud (lu-decompose A allow-zero-pivot))
  79.          (rank (lu-rank (car lud))))
  80.     (if (= rank (* 2 n))        ;ok!
  81.         (lu-backsubstitute (car lud) (cadr lud) B)
  82.         (error "N is too large -- PROCEDURE->RAT-FUNC" n))))))
  83.  
  84.  
  85. ;;; This assembles the matrix of coefficients and the drive vector.
  86. ;;;  The unknowns are in the order a[0], ... a[n-1], b[0], ... b[n-1].
  87. ;;; The rows are produced by evaluating f at 2n+1 test points.
  88.  
  89. (define (assemble-matrix-equation f n)
  90.   (let ((equations (map (lambda (i) (eqn i f n)) (test-points (* n 2)))))
  91.     (list (matrix-by-rows (map cadr equations))
  92.       (list->vector (map car equations)))))
  93.  
  94.  
  95. ;;; Test points are generated along a line parallel to the imaginary axis
  96. ;;; and spaced out to fit within a range near the unit circle.
  97.  
  98. (define (test-points num)
  99.   (map (lambda (imag) (make-rectangular test-point-real imag))
  100.        (let ((space (/ test-point-imag (- num 1))))
  101.      (let loop ((n num))
  102.        (cond ((= n 0) '())
  103.          ((odd? n)
  104.           (cons (* n space) (loop (- n 1))))
  105.          (else
  106.           (cons (* -1 n space) (loop (- n 1)))))))))
  107.  
  108. ;;; The following magic numbers are arbitrary
  109.  
  110. (define test-point-real (/ 1 10))
  111. (define test-point-imag 1)
  112.  
  113. ;;; Returns a representation of an equation:
  114. ;;;  x^n*f(x) = 
  115. ;;;     a[0] + ... + x^(n-1)a[n-1] - f(x)*b[0] - ... - x^(n-1)*f(x)*b[n-1]
  116. ;;; as
  117. ;;; ( x^n*f(x)  ( 1 ... x^(n-1) -f(x) ... - x^(n-1)*f(x) ) )
  118.  
  119. (define (eqn x f n)
  120.   ((with-reasonable-object x
  121.      (lambda (x)
  122.        (values x (f x))))
  123.    (lambda (x fx)
  124.      (let lp ((i 0) (xp 1) (acoeffs '()) (bcoeffs '()))
  125.        (if (= i n)
  126.        (list (* fx xp)               ;rhs
  127.          (append (reverse acoeffs)     ;lhs
  128.              (reverse bcoeffs)))
  129.        (lp (1+ i)
  130.            (* xp x)
  131.            (cons xp acoeffs)
  132.            (cons (* -1 xp fx) bcoeffs)))))))
  133.  
  134. (define ((values . list) receiver)
  135.   (apply receiver list))
  136.  
  137.  
  138. (define (rat-func-compose rf1 rf2)
  139.   (let ((rf1 (check-rat-func rf1))
  140.     (rf2 (check-rat-func rf2)))
  141.     (let ((rf1 (pad-rat-func rf1))
  142.       (rf2 (pad-rat-func rf2)))
  143.       (define (make-list-poly-powers poly)
  144.     (let loop ((power '(1))
  145.            (index 1))
  146.       (if (= index (length (car rf1)))
  147.           (list power)
  148.           (cons power
  149.             (loop (mul-polys power poly) (1+ index))))))
  150.       (let ((numer-power-list (make-list-poly-powers (car rf2)))
  151.         (denom-power-list (make-list-poly-powers (cadr rf2))))
  152.     (define (compose-poly&rat-func poly)
  153.       (apply add-polys
  154.          (map mul-polys
  155.               (reverse (map scale-poly poly numer-power-list))
  156.               denom-power-list)))
  157.     (list (compose-poly&rat-func (car rf1))
  158.           (compose-poly&rat-func (cadr rf1)))))))
  159.  
  160.  
  161. (define (pad-rat-func rf)
  162.   (define (pad-coeff-list coeffs full-length)
  163.     (let loop ((c (reverse coeffs))
  164.            (extra (- full-length (length coeffs))))
  165.       (if (= extra 0)
  166.       (reverse c)
  167.       (loop (cons 0 c)
  168.         (-1+ extra)))))
  169.   (let ((m (max (length (car rf))
  170.         (length (cadr rf)))))
  171.     (map (lambda (coeffs) (pad-coeff-list coeffs m))
  172.      rf)))
  173.        
  174.  
  175. (define (check-rat-func r)
  176.   (if (and (pair? r)
  177.        (pair? (car r))
  178.        (pair? (cdr r))
  179.        (pair? (cadr r)))
  180.       r
  181.       (error "badly formed rational function --" r)))
  182.